Spiral order by disorder and lattice nematic order in a frustrated Heisenberg 
antiferromagnet on the honeycomb lattice 
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Motivated by recent experiments on Bi3Mn40i2(N03), we study a frustrated J1-J2 Heisenberg 
model on the two dimensional (2D) honeycomb lattice. The classical J1-J2 Heisenberg model on 
the 2D honeycomb lattice exhibits Neel order for J2 < Ji/6. For J2 > Ji/6, it has a family of 
degenerate incommensurate spin spiral ground states where the spiral wave vector can point in any 
direction. Spin wave fluctuations at leading order lift this accidental degeneracy in favor of specific 
wave vectors, leading to spiral order by disorder. For spin S = 1/2, quantum fiuctuations are, 
however, likely to be strong enough to melt the spiral order parameter over a wide range of J2/J1. 
Over a part of this range, we argue that the resulting state is a valence bond solid (VBS) with 
staggered dimer order - this VBS is a lattice nematic which breaks lattice rotational symmetry. Our 
arguments are supported by comparing the spin wave energy with the energy of the VBS obtained 
using a bond operator formalism. Turning to the effect of thermal fiuctuations on the spiral ordered 
state, any nonzero temperature destroys the magnetic order, but the discrete rotational symmetry 
of the lattice remains broken resulting in a thermal analogue of the nematic VBS. We present 
arguments, supported by classical Monte Carlo simulations, that this nematic transforms into the 
high temperature paramagnet via a thermal phase transition which is in the universality class of the 
classical 3-state Potts (clock) model in 2D. We discuss the relevance of our results for honeycomb 
magnets, such as Bi3M40i2(N03) (with M=Mn,V,Cr), and bilayer triangular lattice magnets. 



I. INTRODUCTION 

Frustrated quantum magnets support a variety of re- 
markable ground states which emerge as a result of quan- 
tum fluctuations within a large set of classically degen- 
erate configurations.^ Such ground states include valence 
bond solids, magnetic analogues of supersolids, and quan- 
tum spin liquids with various kinds of topological or- 
der. While Neel order is common in bipartite lattices, 
the presence of further neighbor interactions can frus- 
trate this order and lead to interesting quantum ground 
states. This has been extensively studied on the square 
lattice^i"— , and, for S — 1/2, there is an indication of 
a non-magnetic ground state (for 0.45 % J2I J\ ^ 0.6) 
sandwiched between two coUinear magnetically ordered 
ground states. In this paper, we study the J1-J2 Heisen- 
berg model on the honeycomb lattice as the simplest 
model Hamiltonian which incorporates frustration effects 
in this lattice geometry. The Hamiltonian for this model 
is 



7? = Ji ^ S, 



J2 E 



(1) 



where {ij) denotes nearest neighbor pairs of sites, {{ij)) 
denotes next neighbor pairs of sites, and we set Ji , J2 > 0. 

A summary of the results contained in this paper is as 
follows. We find that the classical (5 = 00) model has a 
Neel ordered ground state for J2 < Ji/6. For J2 > Ji/6, 
this gives way to a one parameter family of classically de- 
generate coplanar spin spiral ground states. At 0{1/S), 
quantum fiuctuations within this classical manifold pick 



specific spiral wavevectors, leading to spiral order by dis- 
order. For spin S = 1/2, quantum fluctuations at T = 
are likely to be strong enough to wipe out the spiral order 
parameter over a wide range of J2 / Ji ■ Over a significant 
part of this range of J2/J1, we argue that the spiral or- 
der for spin 5 = 1/2 melts into a valence bond solid 
with staggered dimer order - this state has a spin gap 
and preserves translational symmetry but breaks lattice 
rotational symmetry leading to a 'lattice nematic'. Turn- 
ing to physics at nonzero temperature, spin-spin corre- 
lations decay exponentially at any temperature, but we 
show that the nematic order survives - this nematic trans- 
forms into the symmetric high temperature paramagnet 
via a thermal phase transition which is in the universal- 
ity class of the classical 3-state Potts (clock) model in 
two dimensions. Some of the results on the classical de- 
generacy and spin wave fluctuations have been discussed 
earlier jSiiS, but are included for completeness and clarity. 
We also discuss the connection of our work to previous 
work on this model^rJi and related modelsi^'^'^ 

Before we get into the detailed analysis of the 
above model, we briefly discuss possible materials 
which might realize the physics discussed in this paper. 
Bi3Mn40i2(N03) appears to be an example of a honey- 
comb lattice quantum magnetjlS Since Mn forms MnOg 
octahedral units, and there is strong Hund's coupling, the 
Mn^+ ions behave as S — 3/2 spins. Despite the bipar- 
tite nature of the lattice, and a large antiferromagnetic 
Curie- Weiss constant Qcw ~ —257K, this system shows 
no magnetic order down to T = OAKi^ It has been 
suggested that this arises from frustration due to fur- 
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ther neighbor interactions.'^® Neutron scattering studies 
would be valuable to clarify whether such next neighbor 
couplings are present and whether they place this sys- 
tem in the regime of fluctuating spiral order, leading to 
interesting spin liquid behavior over a wide range of tem- 
peratures, or if there is a nematic transition with its spe- 
cific heat signature being obscured by background lattice 
contributions. Variants of this system, where Mn'*+ is re- 
placed by (with 5=1/2) or by Cv'^+ (with S = 1) 
would also be interesting to study, with the V'*"'" material 
being a possible candidate for observing the dimer solid 
discussed in this paper. 

Among other honeycomb materials, InCu2/3Vi/3 03 
has spin- 1/2 Cu^^ ions nominally forming a honeycomb 
lattice'^ with the nonmagnetic V^"*" ions lying at the cen- 
ter of the honeycomb hexagons. However, this system 
appears to have strong structural disorder since V^"*" and 
Cu^"*" do not order perfectly in this fashion. Ingredi- 
ents for other such honeycomb spin systems could be, 
for instance, Cu^"*" ions within CuOs units arranged on 
the honeycomb lattice. In such a trigonal bipyramidal 
crystal field environment of the oxygens, the copper ion 
would then have a single hole, which is located in the 
d3z2_j.2 orbital. If the resulting S' = 1/2 moments have 
significant next neighbor interactions, they might also be 
candidates to explore the physics discussed here. 

Our study is also relevant to bilayer triangular antifer- 
romagnets where the triangular layers have an AB stack- 
ing, with antiferromagnetic exchange couplings present 
between neighboring sites within each layer (J2) as well 
as between neighboring sites across the two layers (Ji). 
In this case, each layer acts as one sublattice of the 
honeycomb antiferromagnet. Such a structure occurs in 
LuCuGa04 which has copper/gallium ions arranged ran- 
domly in a bilayer triangular lattice leading to a strongly 
disordered spin liquid^ A variant such as HfCu204, if 
it could be synthesized in this structure, might be an 
interesting material to study. 

This paper is organized as follows. We begin, in Sec- 
tion II, with a study of the classical model and its many 
degenerate spiral ground states and follow it up with an 
analysis of spin wave fluctuations and how it selects cer- 
tain spiral ordered ground states from this manifold. We 
argue that spin wave fluctuations are likely to melt the 
order for 5 = 1/2 over a wide range of J2/ Ji- Section III 
contains a bond operator approach to the energetics of 
the nematic (staggered) dimer solid on the honeycomb 
lattice. Section IV describes the effect of thermal fluctu- 
ations on such a nematic state using Landau theory as 
well as by direct Monte Carlo simulations of the classical 
J1-J2 Heisenberg model. Section V contains a discus- 
sion of earlier work on this model and related models 
on other lattices which share some of the features of the 
honeycomb model we have studied. 



II. SPIRAL ORDER FROM QUANTUM 
DISORDER 

A. Degeneracy of coplanar classical ground states 

To calculate the classical ground state energy we be- 
gin by assuming coplanar spiral order on the lattice and 
parameterizing the spins on the two sublattices as 

Si(r) = 5 [cos(Q • r)z + sin(Q • r)x] (2) 
S2(r) = -S'[cos(Q •r + 0)z + sin(Q •r-|-(/))x] (3) 

where Q is the spiral wavevector, r denotes sites on the 
triangular lattice basis, and </) -f- tt is the angle between 
spins on the different sublattices at the same site r. This 
notation is chosen so that the Neel state corresponds to 
Q = (0, 0) and i/i = 0, with spins aligned along ±z. 
The classical ground state energy per spin is given by 

E J 

= Y-i^os(j)+cos{(t)-Qb)+cos{(f)-Qa-Qb)] 

+ J2S^ [cos Qa + COS Qb + cos{Qa + Qb)] , (4) 

where a = x, and b = —x/2 + yy/3/2, are unit vectors 
depicted in Fig. ([1]). Minimizing this classical energy, we 
find that the minimum energy solution for J2/J1 < 1/6 
corresponds to Q* = (0, 0), 0* = 0, so that the Neel state 
is stable for this range of frustration. 

For J2/J1 > 1/6, the minimum energy solutions corre- 
spond to Q* satisfying the relation 



cos Ql + cos Ql + cos(g* +Q*b)^ 




while 0* is determined completely by 

sin0* = 2:^(sing^-Hsin(Q: + Q,*)), (6) 
COS0* = 2^(H-cosg^ +cos(Q: + Qn). (7) 

It is clear that the spiral wavevector is not uniquely 
fixed by the above relations, as has been noted earlier i^iiS 
As shown in FiglU the set of classically degenerate so- 
lutions to Eq.lO (which we label Q*) form a closed 
contour^ around Q=(0,0) for 1/6 < J2/J1 < 1/2. For 
^2/^/1 > 1/2, it forms closed contours around (Qa, Qb) = 
±(27r/3, 27r/3). This regime does not appear to have been 
investigated in earlier work. In the limit J2/J1— >oo, 
where the two triangular sublattices of the honeycomb 
lattice approximately decouple, Q* ±(27r/3, 27r/3) 
which is the ordering wavevector of the 120° state on 
the triangular lattice. We focus next, therefore, on how 
quantum or thermal fluctuations select specific spin spi- 
rals from the manifold of classically degenerate coplanar 
spirals discussed above. 



FIG. 1: (Color online) Left panel: Real space basis vectors 
for the honeycomb lattice. Right panel; Momentum space 
picture depicting the manifold of classically degenerate spiral 
wavevectors for J2/J1 =0.3 (red, thin solid), J2/J1 = 0.5 (pur- 
ple, dash-dotted), and J2/ Ji=0-7 (green, dashed). Also indi- 
cated by purple (solid) dots are the six distinct spiral wavevec- 
tors lying on this manifold which are favored by quantum 
fluctuations. Black (thick solid) hexagon indicates the first 
Brillouin zone of the lattice. 



B. Spin wave fluctuations 

We calculate leading quantum corrections using 
Holstein-PrimakofF (HP) spin wave theory. We begin by 
defining new spin operators S via 




'cos6'f(r) —s'm9g{r 

1 
sin0^(r) cos0^(r) 




(8) 



where £ = 1,2 labels the sublattice, 9i{r) = Q • r, and 
62{r) = Q-r + (j). Reexpressing the Hamiltonian in terms 
of these new spin operators and rewriting these spin op- 
erators in terms of HP bosons, we arrive at the following 
Hamiltonian which includes the leading spin wave cor- 
rection to the classical ground state energy. 



+ 2S 



k>0 



pi 



2Al 



Here &t = (^{(k) 4(k) &i(-k) b^{-k)), Ek>o indicates 
that the sum runs over half the first Brillouin zone (so 
that k and — k arc not both included), and the Hamilto- 
nian matrix Mk takes the form 



Bl Ak Dl Ck 

Ck £'k ^k ^k 

\DT Ck BT AJ 



(10) 



with explicit expressions for A\^-D\^ given in Appendix 
A. Diagonalizing this problem via a generalized Bogoli- 
ubov transformation, we obtain the spin wave corrected 
ground energy as 



Ea + 2SJ2 [^-(k) + A+(k) - 2Ak] 

k>0 



(11) 



FIG. 2: Plot of the spin wave correction to the energy per site 
AE — {Equ — Eel) /N (in units of Ji) as a function of J2/J1. 



The eigenvalues A±(k) are given by 



A±(k) = v/ak ± I3k 



where 



Qfk = A 



/3k 



(12) 



(13) 



MAkB^^-CM^ + iDkB^-B^D^y. (14) 



For J2 = 0, it is known from quantum Monte Carlo 
simulations that this model has long range Neel order,— 
We have checked that the Neel state energy for S* = 1/2 
is, for J2 = 0, in good agreement with recent quantum 
Monte Carlo simulations in the valence bond basis. 

The quantum correction to the classical ground state 
energy is responsible for selecting a unique quantum 
ground state from the manifold of classically degenerate 
ground states. Minimizing this energy correction over 
the classical ground state manifold, Q*, we find the fol- 
lowing results for the spiral wavevector Q**, which is 
selected by quantum fluctuations, with the resulting 
being determined by Eqns. (|6l7l) . 

For 1/6 < J2/J1 < 1/2: The ground state is a spiral 



(9) state 5*1, with 



■ 16 J| 



a' 



= cos""(TT-bT - -) (15) 

QT = (16) 

While the above relations specify a single spiral state, 
there are a total of six symmetry related spirals, the other 
five being obtained by 27r/6 rotations of the above Q**. 

For 1/2 < J2/J1 < 00: The ground state is a different 
spiral state 52, with 

Ql* = TT — cos" 



'(— + -) 



1. 



2 cos 



(17) 
(18) 



^4J2 ^ 2> 

There are six symmetry related ^2 spirals, the other five 
being obtained by 27r/6 rotations of the above Q**. The 
spin wave correction to the ground state energy is shown 
in Fig. dH). 
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C. Spiral order parameter 'melting' 

Spin wave fluctuations will tend to renormalize the spi- 
ral order, and may render the spiral states unstable. We 
have checked that the leading spin wave correction to the 
spiral order parameter, given by 

]^ w) = -wj^ E(^l (k)&,(k)), (19) 

diverges as log(A'^) since the spin wave energy vanishes on 
the entire classical manifold of degenerate spiral wavevec- 
tors. This suggests that the spiral order will disappear 
for any value of S. However, while such line zeroes of 
the dispersion ('Bose surfaces') are mandated by con- 
servation laws in the compressible phase of certain ring- 
exchange models, it is not generic in this model, and the 
spin wave energy is expected to only have gapless points 
in momentum space corresponding to those wavevectors 
which are selected by quantum fluctuations. We expect 
spin wave interactions, not included at this stage, to gap 
out all other wavevectors and stabilize the spiral order 
for large enough S. This requires a higher order spin 
wave calculation (in powers of 1/S) which is beyond the 
scope of this paper. At this stage, we restrict ourselves 
to noting that an exact diagonalization studj*iS of the 
spin-1/2 model did not find any evidence of a tendency 
towards magnetic ordering over a wide range of J2/J1 
where the classical analysis predicts spiral order. Based 
on this, we expect that while spiral order might be sta- 
bilized at large S from spin wave interactions, this order 
is likely to 'melt' for small spin values, leading to other 
competing states. 



III. FLUCTUATION INDUCED 'LATTICE 
NEMATIC ORDER 

We have seen that quantum fluctuations in a spin wave 
expansion will tend to strongly suppress and, for small 
spin values, perhaps disrupt the spiral order. In accor- 
dance with the Mermin- Wagner theorem, thermal fluc- 
tuations are similarly expected to melt the spiral order 
for any nonzero temperature. While such quantum and 
thermal fluctuations may restore spin rotational symme- 
try with exponentially decaying spin correlations, there 
could be persisting broken symmetries in bilinears of the 
spin operator (which are obtained, for instance, by taking 
dot products or cross products of the single spin opera- 
tors). We begin by listing such bilinears in order to see 
which of them could possibly survive the effect of fluctu- 
ations that destroy magnetic long range order. 

In the ordered spiral state, ignoring spin wave correc- 
tions to the correlation functions, we find, 

Si(r)xSi(r+R) = S2(r)xS2(r-|-R)=52sin(Q-R)t/(20) 
Si(r)xS2(r±R)= -S'2sin(Q • Ii±(t>)y. (21) 



Such bilinears therefore preserve lattice translational 
symmetry but break the rotational invariance of the lat- 
tice. Since solutions (Q, </>) and (— Q, — </>) are related by 
a global spin rotation, and spin correlations are short- 
ranged at nonzero temperature, such 'vector chiralities' 
are also expected to have exponentially decaying corre- 
lations at nonzero temperature. By contrast, spin corre- 
lations such as 

Si(r)-Si(r+R) = S2(r)-S2(r-|-R)=S'2cos(Q-R) (22) 
Si(r)-S2(r±R)= -52cos(Q-R±0) (23) 

are invariant under global spin rotations. These correla- 
tions are clearly invariant under lattice translations, but 
they break lattice rotational symmetry. Such a discrete 
broken symmetry may survive even after fluctuations ren- 
der the spiral state unstable. 

Let us focus on nearest neighbor bonds and write out 
the above spin correlations which are simply proportional 
to the bond energies. We find, for the three bonds around 
a site on sublattice-1, 

Si(r)-S2(r) = cos (24) 
Si(r)-S2(r-6) = -S^ cos{Qi, - (b) (25) 
Si(r)-S2(r-a-6) = -S^ cos{Qa + Qt - ^) (26) 

Computing these bond energies in the Si spiral ground 
states selected by quantum fluctuations, we find that two 
of the three bond energies are equal while the third takes 
on a different value, so that the C3 rotational symmetry 
about a lattice site is broken in the Si state. This is the 
three-fold symmetry that we expect may still be broken 
even if spin rotational symmetry is restored by quantum 
or thermal fluctuations. The S2 state also breaks the 
three-fold lattice rotational symmetry, as seen from the 
corresponding spiral ordering wavevectors. Fluctuations 
about the spiral states could thus lead to a 'lattice ne- 
matic' state, which is invariant under lattice translations 
but not lattice rotations. Below, we discuss a quantum 
nematic valence bond solid (VBS) state as a candidate 
ground state for S = 1/2, as well a classical nematic state 
induced by thermal fluctuations for any spin value. 

A. Nematic valence bond solid 

Motivated by the above discussion, we consider the 
simplest candidate for a lattice nematic ground state 
for S' = 1/2 spins, which corresponds to forming a ne- 
matic Valence Bond Solid (VBS) which consists of singlet 
dimers on the honeycomb lattice as shown in Figl3] Such 
a state has been proposed earlier over a small window of 
J2/J1 on the basis of a small system exact diagonaliza- 
tion study An amusing way to view this VBS state, as 
shown in Fig. ([3]), is to think of it as arising from cou- 
pling together frustrated spin 5 = 1/2 J1-J2 chains. If 
we imagine the interchain couplings being tuned to zero, 
this would lead to decoupled Majumdar- Ghosh chains 
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constraint 



FIG. 3: Sketch of the valence bond solid state with staggered 
dimer order which breaks the honeycomb lattice rotational 
symmetry but preserves spin rotational and lattice transla- 
tional symmetries, a, b denote basis (unit) vectors of the tri- 
angular lattice formed by the dimer bonds, di.2 and Ai,2 
indicate Hartree-Fock-Bogoliubov mean field parameters in 
the bond operator mean field theory (see text for details). 
Dashed line indicates the set of spins which might be viewed 
as forming one dimensional dimerized chains which are cou- 
pled in the transverse direction. 



E 



t"! t 
r a r,a 



(31) 



Here a is summed over x, y, z, ta being the three triplon 
operators on the bond at r. (Repeated Greek indices 
henceforth denote summation over x,y,z.) To simplify 
the calculation, we satisfy this constraint on average, 
rather than locally. We assume the singlet operators to 
be condensed, allowing us to replace the operator Si with 
a number s. The excitations are the triplet operators, 
and the terms of the Hamiltonian may now be organized 
in order of the number of triplet operators. The Hamil- 
tonian in momentum space, keeping only terms up to 
quadratic order, is 



H 



[2] . 
BO' 



k>0 



k>0 



Gk Fk 




■ i7(k) " 


_F* Gk_ 




_tt(-k)_ 



(32) 



which are known to possess dimer order with a spin gap; 
in particular, the dimerized state is the exact ground 
state of the single chain at J2 = 0.5Ji. The honeycomb 
lattice VBS might then be thought of as arising from the 
decoupled chain limit upon incorporating interchain cou- 
plings while leaving the singlet gap intact. The choice 
of which direction these chains run along is completely 
arbitrary in the honeycomb limit, so that there are three 
degenerate ground states that break the G3 lattice rota- 
tional symmetry. We note that Heisenberg models with 
multispin interactions have been proposed for which this 
VBS state is the exact ground statei^ 

In a state where such singlets are forced to occur on 
the indicated bonds, the only excitations correspond to 
breaking these singlets to form triplet excitations which 
can then form a 'triplon' band. To analyze the energetics 
and stability of such a state, we therefore use the bond 
operator formalism proposed in Ref. [2^. The details of 
the calculation are presented in Appendix B. 

Instead of working in the nai've basis of S — 1/2 op- 
erators on every site, we switch to a basis of singlet and 
triplet bosonic operators defined as 



5t|0) = 




n)- 


Ut)) 


(27) 


4lo) = 


^< 


n)- 


\m 


(28) 


4lo) - 


T.^ 


n) + 


\m 


(29) 


4lo) - 


T.^ 


n) + 


ut)) 


(30) 



on the dark (dimer) bonds in Fig|31 together with a local 



where is a chemical potential that has been introduced 
to satisfy the constraint in Eq. (|31[) . The matrix entries 
are given by 

Ji 

Gk = ^ - ^•^i(ek + e-k) + ^J2(?7k + ?/-k) 
Fk = -^^i(ek + e-k) + J2^(?7k + '7-k) 
where 



77k = 2[cOs(fca) + COs(fcfc) -I- COs(fca -I- fcfe)] 

Diagonalizing the Hamiltonian by a bosonic Bogoliubov 
transformation gives the dispersion of the eigenmodes to 
be 

Ek - ^Gl-\Fk\^ (33) 

The energy of this state is plotted as the green dashed 
line in Fig (g]). 

While this quadratic theory gives a consistent picture 
of our lattice nematic state, higher order terms may lower 
its energy significantly. We proceed to take these into 
account by means of a self consistent Hartree-Fock ap- 
proach. This approach has been compared recently, for a 
star lattice Heisenberg antiferromagnet, with Gutzwiller 
projected variational wavefunctions^^ and exact diago- 
nalization studies^^ and shown to provide a good descrip- 
tion of the energetics of valence bond solid states on the 
star lattice.— We begin by noting that the terms of cu- 
bic order in the triplet operators do not contribute, since 
we work with the assumption that the triplon operators 
themselves are not condensed. The quartic part of the 
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J2/J1, we expect competing states might emerge which 
are descendants of spin hquid states on the triangular 
lattice^^ - this needs further investigation. 

We note that this bond operator formahsm does not 
take into account the fluctuations of the singlets them- 
selves; the kinetic energy lowering from such resonat- 
ing singlet valence bonds might possibly favor i/S x 1/3 
plaquette dimer order which also breaks lattice transla- 
tion symmetryii^S - such states can be accessed within a 
Schwinger boson formalism and might be relevant in the 
vicinity of the point where the Neel order is lost. 



FIG. 4: (Color online) Ground state energy (in units of Ji) 
as a function of J2/J1. The red (solid) line is the energy of the 
spiral state including leading order spin wave corrections, the 
green (dashed) line is the nematic VBS energy up to quadratic 
order in triplon operators, and the blue (dash-dotted) line 
indicates nematic VBS energy up to quartic order in triplon 
operators. 



Hamiltonian is given by 

k,k',q 

tl{k + q)mk' + q)tl{k')ts{k) (34) 

where e^^g-y is the permutation symbol. Guided by the 
symmetry of the nematic phase, we postulate the follow- 
ing real-space order parameters: 



2 (^r,7^r+(Si ,7) 

2 (^r,7^r+(52,7) 
1 

— (ir,7^r+(5i ,7/ 



di 
Ai 

A2 = ■:j(ir,7^r+(52,7/ 



(35) 
(36) 
(37) 
(38) 



where = ±a, and 82 — ±6,±(a-|-6). These are defined 
on bonds as shown in Fig. ([3]). The above are the only 
operators that couple to at quadratic level. 

We calculate these order parameters self-consistently, 
and thereby obtain the energy of the nematic VBS, hav- 
ing accounted for quartic terms. This is plotted in 
Fig. @ as the blue (dot-dashed) line. At quadratic level, 
the nematic state is energetically favourable over the spi- 
ral over a small window near J2 ~ 0.35Ji. Although 
we have not considered interactions between spin wave 
modes, it is nevertheless, it is encouraging that the quar- 
tic level energy in the bond operator formalism is lower 
than the spin wave energy for J2 ^ 0.25Ji, except for a 
small window around J2 — 0.5Ji. Since the spiral order 
is anyway likely to be suppressed by fluctuations, our re- 
sults are quite suggestive of such nematic VBS order be- 
ing present over a wide window of frustration. At large 



B. Thermal fluctuations: Landau theory 

The spiral states Si and S2, obtained from including 
spin fluctuations at large S, can only be stable at zero 
temperature. At any non-zero temperature, since our 
system is two dimensional, spin rotational symmetry will 
be immediately restored. As earlier discussed, the sim- 
plest ordering would involve nearest-neighbor bilinears 
of the spin operators, which may break lattice rotational 
symmetry. Motivated by earlier work on such 'lattice 
nematics' and quantum dimer models, we define a local 
complex order parameter 

^(r) = (Si(r).S2(r))+c^(Si(r).S2(r-6)) 

+ cj2(Si(r)-S2(r-a-6)) (39) 

on sites of sublattice 1, where w — exp(z27r/3). Since 
two of the bond energies are equal in the ground state, 
'(/'(r) ~ {1,0;, o;^} in the three ground states, so that 
'4>^{v) ~ 1. This order parameter is invariant under 
translations even in the spiral state. Under an anticlock- 
wise 27r/3 rotation about a site on sublattice 1, we have 
ijj — >■ ioif). Finally, reflections about axes running along 
the bonds can be shown to lead to V — > V'*- Based on 
these symmetry considerations, the finite temperature 
classical lattice nematic is expected to be described by 
a Landau free energy functional of the form 



F = J d^r [m\^P{r)\^ + u\l|Jir)\^ + X\Vij{r)\'^ 

+ w{^P^r) + tl;*^(r))] (40) 

which is equivalent to a 3-state clock model (or equiv- 
alently, the 3-state Potts model) if we assume that am- 
plitude is fixed. We therefore expect the classical model 
(as well as possibly the 5=1/2 model) to exhibit, upon 
warming up from T = 0, a finite temperature transi- 
tion from a lattice nematic into an ordinary paramagnet, 
with this transition being in the universality class of the 
3-state Potts model in 2D. (Of course, these arguments 
do not rule out the possibility of a first-order transition.) 
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FIG. 5: Binder cumulant of the order parameter (from 
Eq. I4ip plotted as a function of temperature for various sys- 
tem sizes. The crossing point of the curves at T / J\S^ « 
0.1515(5) indicates a continuous nematic to paramagnet phase 
transition at J2/J1 =0.8. Inset shows the nematic transition 
temperature, Tc, as a function of J2/J1. Lines are guides to 
the eye. 



C. Thermal fluctuations: Monte Carlo study 

We have carried out Monte Carlo simulations of the 
classical J1-J2 Heisenberg model in order to numerically 
explore the nematic-paramagnet phase transition. Using 
a combination of single-spin Metropolis moves and energy 
conserving (over-relaxed) moves, we have computed the 
Binder cumulant of the order parameter. 



B = 1 



1*1 



3(1*1 



(41) 



where * — J^r'^i^) ^ complex scalar, and the suscep- 
tibility. 
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FIG. 6: Nematic susceptibihty (from Eq. [42]), in units of 
1/Ji, as a function of temperature for various system sizes 
at J2/J1 = 0.8. Lines are guides to the eye. Inset shows the 
specific heat peak at the transition for L = 72. 



reasonably consistent with a 3-state Potts model tran- 
sition for which the exact exponents are known^^ to be 
/3 — 1/9,7 — 13/9, and v = 5/6; these imply ^jv ^ 1.733, 
fijv K, 0.133, and Xjv — 1.2. The critical Binder cumu- 
lant B{Tc) also seems consistent with earlier numerical 
work on the Potts model<^ We have also computed the 
specific heat of this model, and, as seen in the inset of 
Figini it shows a clear peak at the transition point lo- 
cated from the Binder cumulant calculation. For reasons 
we do not completely understand, the specific heat does 
not exhibit clear finite size scaling over the system sizes 
explored. It is possible that the thermally fiuctuating 
spin wave modes and their interaction with the nematic 
order parameter may make it difficult to extract the fi- 
nite size scaling of the specific heat singularity for system 
sizes we have studied. We are carrying out careful nu- 
merical studies of this model on larger system sizes in 
order to understand this issue. 



(42) 



for the classical Heisenberg model for various J2/Ji- 
Fig. [5] shows the Binder cumulant as a function of tem- 
perature obtained on various system sizes {N — with 
L = 36-72) for J2/J1 = 0.8 by averaging over 10^- 
10^ configurations. These exhibit a crossing point at 
Tcj J\S^ « 0.1515(5) indicating a continuous thermal 
phase transition, with B{Tc) « 0.63. 

In addition, as seen from Fig.[6l the peak of the suscep- 
tibility (at J2/J1 =0.8) increases with system size. Based 
on the finite size scaling of this peak height, x ~ L''^'^, 
we find 7/j/ « 1.68(8). The order parameter (|*|) at 
scales with system size as L~^^'^, with j3/v 0.14(2). 
Finally, the shift in the susceptibility peak with system 
size is expected to scale as AT-^ ^ L~^/'^ from which we 
find l/v « 1.25(9). These results for the exponents are 



IV. RELATION TO PREVIOUS WORK 

Earlier investigations of the honeycomb lattice model 
have focused on the spin wave selection of various spi- 
ral statesi^iiS Our results are in line with these stud- 
ies - it appears that specific spin spirals are selected at 
0(1/5) in a spin wave calculation, but the resulting order 
is likely to 'melt' for 5' = 1/2 over a wide range of J2/J1. 
An exact diagonalization study of the spin S* = 1/2 
model has suggested that nematic order with breaking of 
C3 rotational symmetry could appear in the vicinity of 
J2 =0.4-0.5 Ji, and this order has also been guessed from 
a study of Berry phase effects in a nonlinear sigma model 
formulation^^ Our bond operator calculations lend sup- 
port to this claim, and also suggest that this nematic 
dimer order may persist over a wide range of J2/J1. 
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The idea that isotropic Heisenberg models may have 
such nematic orders at finite temperature is well known 
from early work on the square lattice J1-J2 model.— 
For J2 < Ji/2, the classical {S = 00) ground state of 
this model on the square lattice is Neel ordered. For 
J2 > Ji/2, by contrast, there is a large set of classi- 
cally degenerate ground states in which the two sublat- 
tices are individually perfectly Neel ordered with an arbi- 
trary relative angle between the two sublattices. Within 
this classical manifold, quantum fluctuations at 0{1/S) 
select coUinear ground states^i with ordering wavevec- 
tors Q ~ (tt, 0) or (0,7r). At any nonzero temperature, 
this model exhibits exponentially decaying spin correla- 
tions, consistent with the Mermin- Wagner theorem, but 
the broken lattice rotational symmetry associated with 
these coUinear ground states survives at low tempera- 
ture. Upon further heating, this 'lattice nematic', which 
breaks the C4 rotational symmetry of the square lattice 
down to C2, converts into the high temperature param- 
agnetic phase via an Ising transitional^ Despite a large 
number of numerical studies however, the ground 
state phase diagram of this square lattice spin-1/2 model 
appears to not to be satisfactorily understood. 



The relation between spiral magnetic states and ne- 
matic orders has also been explored in the context of the 
J1-J3 model on the square lattice^2S In this case, there is a 
Neel to spiral transition for J3 > Ji/4, which is a Lifshitz 
transition similar to the case we have studied. Melting 
this spiral thermally leads to an Ising nematic similar to 
the square lattice J1-J2 model. The main differences of 
our model with this case are: (i) The spiral wavevector 
is unique (modulo reflections) in the classical square lat- 
tice J1-J3 model unlike the line degeneracy we encounter 
on the honeycomb lattice; (ii) the nematic-paramagnet 
transition in the square lattice J1-J3 model is an Ising 
transition; (iii) unlike on the honeycomb lattice, there 
is no simple quantum analogue of the classical nematic 
in the square lattice model. It may be more useful to 
consider possible analogies of the honeycomb model with 
the J1-J2-J3 model on the square lattice which has been 
studied in recent worki^ 



V. SUMMARY 



We have studied the honeycomb lattice J1-J2 Heisen- 
berg model. We have seen that the classical model sup- 
ports a one-parameter family of degenerate spin spiral 
states, of which specific spin spirals are selected out by 
quantum fluctuations. For general spin values, we expect 
the spiral order to be strongly suppressed but robust ne- 
matic order to survive. For S = 1/2, spin fluctuations are 
likely strong enough to 'melt' the spiral order leading to 
a spin gapped nematic dimer solid as indicated from our 
bond operator calculations. We have shown that the clas- 
sical model, and possibly also the dimer solid, are con- 
nected to the high temperature paramagnetic phase via a 
3-state Potts model transition. Neutron scattering exper- 
iments would be valuable to test for fluctuating spiral or- 
der at finite temperatures - in this regime, the equal time 
structure factor exhibits peaks on the spiral contours 
in Fig. [1] as the system thermally explores the various 
(nearly) degenerate spirals. This may allow a determi- 
nation of the further neighbor couplings which frustrate 
Neel order. Further work is necessary to determine if in- 
teresting gapless spin liquids emerge as candidate ground 
states for this model over some regime of frustration as 
has been recently proposed for other frustrated quantum 
magnetsj22r— The other interesting possibility is the ex- 
istence of gapped spin liquids as have been recently un- 
covered in numerical studies in the insulating state of the 
honeycomb lattice Hubbard model. The model we have 
studied appears to be directly applicable as an effective 
spin Hamiltonian (with J2/J1 ~ 0.1) in the insulating 
phase of the Hubbard model for moderate values of re- 
pulsion. Finally, our results are relevant to honeycomb 
and bilayer triangular magnets; we hope our work stimu- 
lates further experiments on the Bi3M40i2(N03) family 
of materials, and other compounds which might realize 
this model. 



Some features of the honeycomb model, such as a 
highly degenerate set of classical spiral states and the 
resulting spiral selection by fluctuation effects, bear 
similarities with studies on the diamond lattice J1-J2 
modelf^i^^ which were motivated by insulating spinel 
compounds such as MnSc2S4, C03O4, and CoRh204. We 
note, in passing, that the issue of spiral order and its 
connection to lattice nematicity also arises in the con- 
text of itinerant systems. In particular, such spiral melt- 
ing has been proposed as one mechanism^ for the ob- 
served nematic transport^ in Sr3Ru2 07 at intermediate 
magnetic fields, although there are competing theoretical 
proposals'^^ for the observed nematic order. 
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Appendix A: Matrix Elements for the 
Holstein-Primakoff Hamiltonian 

The explicit expressions for the matrix elements of the 
Holstein-Primakoff Hamiltonian in Eq. (|10[) are given by 

Ak = ^[cos(j)+cos{(t)-Qb) + cos{<j)-Qa-Qb)] 

- J2[C0S Qa + cos Qb + COs{Qa + Qb)] 

+ -y[(cosQa + l)cosfcQ + (cosQb + l)cosfcb 

+ (cOs(Qa+Q6) + l)cOs(fca+fc6)] (Al) 

Bk - ^[(cos0-l) + (cos(0-Qb)-l)e-^'=^ 

+ (cos((/) - 0„ - Qb) - l)e-*(*'"+^-^)] (A2) 
Ck = Tk + K (A3) 

Tk = ^[(cos(gj-l)e''^-" + (cos(Qb)-l)e-*"-'' 

+ (cos(Qa + g;,)-l)e^(""+'='')] (A4) 

£'k = ^[(cos0+l) + (cos((/)-Qb) + l)e"'''''' 

+ {cos{cj)-Qa-Qb) + l)e-'^''^+''''>] (A5) 

Appendix B: Triplon Calculation for Lattice 
Nematic State 

We work with a basis of singlet and triplet operators 
that are centred on bonds indicated in Fig. In terms 
of the bond operators defined in the text, the spin oper- 
ator on any site can be written as 

S]{r) = lfm4tr,^ + tUsr) - '^e^f!stlpU,s (Bl) 

Here, the index £ indicates the sublattice. The factor f{£) 
takes the value -1-1 on sublattice 1 and —1 on sublattice 
2. r is summed over sites of the bond-centred triangular 
lattice. 

We now consider the singlets to have condensed, giving 
us a nematic state. Between sites that are connected by 
a bond, we have 

Si(r).S2(r) = -j52 + i^4_^ir,^ (B2) 
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For sites that are not connected by a bond, we have 
S,(r).S,, (r') = / + ir,^) X 

To enforce the constraint on the bond operators, we 
rewrite the Hamiltonian as 

i/-/i^[s'+tLir,a-l] (B4) 
r 

where H is the original spin Hamiltonian in Eq. ([1]). /i is 
now tuned so that the constraint is satisfied on average. 
Keeping terms to quadratic order in the i-operators, this 
Hamiltonian may be rewritten as Eq. p2p , and diagonal- 
ized by a Bogoliubov transformation. For fixed J2/J1, 
we choose the value of s that minimizes energy. 

The term that is quartic in triplon operators is given by 
Eq. (IMl) . This can be decoupled in hopping and pairing 
channels, using the order parameters defined in Eq. (j38l) . 
This modifies the coefficients of the quadratic Hamilto- 
nian of Eq. (|32|) as follows 

= Gk + yd2(ek + e-k) + J2rf2(ek + e-k) 

+ 2J2diCOs{ka) (B5) 

4" = i^k-yA2(ek + e-k)- J2A2(ek + e-k) 

- 2J2Aicos(fcJ (B6) 

In addition, the Hamiltonian acquires a constant contri- 
bution given by 

<5^W ^ ^3j^N{dl - lAap) - 6J2N{dl - \A2\'') 

-3J2N{dl^\A,\^) (B7) 

The Hamiltonian is solved by a Bogoliubov transforma- 
tion. For fixed s, the d and A order parameters are de- 
termined self-consistently, while is tuned to make sure 
the constraint on bond operators is satisfied, s is chosen 
to minimize the ground state energy for every J2I Ji- 
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